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Abstract 

We calculate analytically the largest Lyapunov exponent of the so-called aXY 
Hamiltonian in the high energy regime. This system consists of a d-dimensional 
lattice of classical spins with interactions that decay with distance following a power- 
law, the range being adjustable. In disordered regimes the Lyapunov exponent can 
be easily estimated by means of the "stochastic approach", a theoretical scheme 
based on van Kampen's cumulant expansion. The stochastic approach expresses 
the Lyapunov exponent as a function of a few statistical properties of the Hessian 
matrix of the interaction that can be calculated as suitable microcanonical averages. 
We have verified that there is a very good agreement between theory and numerical 
simulations. 
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1 Introduction 

The largest Lyapunov exponent (LLE) measures the sensitivity to initial con- 
ditions in dynamical systems. More precisely, it gives the asymptotic rate of 
exponential growth of most vectors in tangent space. Thus, the determination 
of the LLE, either numerical or analytical, requires the knowledge of the long- 
time behavior of a system of linear differential equations with time-dependent 
coefficients. 

From a numerical point of view, there are standard and simple algorithms for 
computing the LLE that perform satisfactorily, even for large systems [1]. To 
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the best of our knowledge, the first attempt to evaluate analytically the LLE 
of many-body systems with smooth Hamiltonians was made by Pettini and 
collaborators, less than ten years ago [2,3]. This "geometric" approach begins 
by postulating a mean-field approximation in tangent space. Then, the few 
resulting equations are handled with the techniques of stochastic differential 
equations (cumulant expansion [4]) which, eventually, lead to an analytical 
estimate of the LLE. Despite the success of the geometric method in repro- 
ducing the LLE of some many-particle systems, it must be stressed that, from 
a theoretical point of view, this method is unsatisfactory: it presents draw- 
backs such as ad-hoc approximations, ill-defined ingredients, and free param- 
eters (see the discussion in Ref. [5]). However, it is possible to construct a 
first-principles theory following the steps of Pettini et al., but avoiding the 
defects of the geometric method. The idea is simple: do not invoke a mean- 
field approximation at the start, just apply the cumulant expansion to the 
full A^-particle problem. This recipe was partially executed by Barnett et al. 
[6] and then complemented by the present authors [5]. The result is a theory 
where all the ingredients are well defined, there are no free parameters, and 
approximations are controllable. Even if several of these approximations are 
really crude, it was shown [7] that the theory is still capable of quantitative 
predictions. 



Reference [7] contains the application of the "stochastic approach" to the 
disordered regimes of the infinite-range XY Hamiltonian (HMF) [8,9]. This 
model consists of a rf-dimensional lattice of classical spins rotating in parallel 
planes, and with a pairwise interaction that depends on the relative orienta- 
tion of spins but not on the distance between them. In Ref. [10] the HMF 
was generalized by making the interaction depend on the distance r accord- 
ing to the power law r~", thus obtaining the family aXY (the case a = 
recovers the HMF, a/d ^ oo corresponds to first neighbors interactions). The 
purpose of this paper is to estimate analytically the LLE of the disordered 
regimes of the aXY family, using the stochastic approach [7] . The restriction 
to these regimes is not an intrinsic limitation of the theory but obeys to the 
general technical difficulty of determining two-time correlation functions. In 
the disordered regimes of the aXY models the relevant correlation functions 
are approximately Gaussian. 



The rest of the paper is organized as follows. For reasons of self-containedness, 
we begin by presenting a very short review of the theory in Sect. 2. The 
systems to be studied (the aXY Hamiltonians) are introduced in Sect. 3, 
where we also work out the predictions of the theory for these particular 
systems. Section 4 contains the critical comparison of theoretical results with 
numerical simulations. We close in Sect. 5 with a summary and some remarks. 
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2 A short review of the theory 



We consider Hamiltonians of the type 
1 ^ 

n = -Y.p' + nQi,---,<iN), (1) 

^ i=i 

where qi and pi are conjugate coordinates, and the interaction potential V is 
assumed to be smalL This is the case of the disordered phases of the aXY 
models where the interaction energy is very small as compared to the kinetic 
energy. Let us denote phase space points by x and tangent vectors by ^. 
Differentiating the Hamilton equations, one obtains the evolution equations 
for tangent vectors: 

e = Mt) C (2) 



For a Hamiltonian of the special form (1), the operator A has the simple 
structure 

AW=(_4) J). (3) 



Here V is the Hessian matrix of the potential V, namely Vij = d'^V / dqidqj . 
Once initial conditions xq and have been specified, the LLE A can be found 
by calculating the limit [1] 

A - limlln|e(t;a;o,eo)r. (4) 



We assume that for any initial condition xq in phase space, the trajectory 
x{t] Xq) is ergodic on its energy shell. This implies that A depends only on 
energy and other system parameters, but not on xq, which can then be chosen 
randomly according to the microcanonical distribution. There will also be no 
dependence on initial tangent vectors, because if is also chosen randomly, it 
will have a non-zero component along the most expanding direction. Moreover, 
assuming weak intermittency, one can write 

A«^ lnn^ln(|e(t;xo,eo)r) , (5) 



brackets meaning microcanonical averages over xq. This is our first approxi- 
mation, to be discussed in Sect. 5. 

By letting xq be a random variable, V(t; Xq) becomes a stochastic process, and 
Eq. (2) can be treated as a stochastic differential equation. However, as we 
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arc interested in the square of the norm of ^, we focus, not on Eq. (2) itself, 
but in the related equation for the evolution of the "density matrix" 



^^m-A^e+^fA^^Me: (6) 



the rightmost identity defining the linear superoperator A. For the purpose of 
the perturbative approximations to be done, the operator A is split into two 
parts A = Aq + Ai{t), where Aq corresponds to the evolution in the absence 
of interactions, and Ai depends on the Hessian of the interaction. Assuming 
that the fluctuations of Ai(t) are small and/or rapid enough, it is possible to 
manipulate Eq. (6) to derive an explicit expression for the evolution of the 
average of 

{^e){t) = exp{tA) ^0^^ , (7) 



where A is a time-independent superoperator given by a cumulant expansion, 
which we truncate at the second order (the lowest non-trivial order). We refer 
the reader to [4,5] for explicit expressions. The truncation of the cumulant ex- 
pansion is our second approximation. The perturbative parameter controlling 
the quality of the truncation is the product of two quantities, the "Kubo num- 
ber" (TTc. The first factor, a, characterizes the amplitude of the fiuctuations 
of Ai{t). The second, Tc, is the correlation time of Ai{t). 

Let -Lmax be the eigenvalue of A which has the largest real part. Taking the 
trace of Eq. (7), one sees that the LLE A is related to the real part of L^ax by 

A = iRe(L^ax) • (8) 



To proceed further one needs the matrix of A in some basis. The crudest 
approximation consists in restricting A to a particular three-dimensional sub- 
space, a choice that is equivalent to a mean-field approximation in tangent 
space [5]. The corresponding 3x3 matrix for A is: 

/ 

A= 2aVii) -2a^TP -2/x | , (9) 
V-/. + 2aVi2) 

with the definitions 



/.= 1tV(V), a' ^ j^Trmf) , (10) 
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oo 

rf = / dTT'fir) , fir) = ^Tr(5V(0)5V(r)) (11) 



(see [5,7] for details). In the mean-field approximation the LLE is expressed in 
terms of the set of four parameters and a^T^''~^^\ k = 0,1, 2. The parameters 
fj, and a are, respectively, the mean and variance of the stochastic process V(t). 
The characteristic time t^^^ = Tc is naturally interpreted as the correlation 
time of the process V(t), /(r) being the relevant correlation function. 



3 The aXY Hamiltonians 

In this section we apply the perturbative/mean-field theory of Sect. 2 to the 
aXY Hamiltonians 

^f = i^pj + if;i^i2?W^,K + v, (12) 



where N — N{a, d) = J2ij ^ij^ 1^-, and is the distance on the periodic lattice 
[10] . This model represents a lattice of classical spins with the interaction range 
controlled by the parameter a. Each spin rotates in a plane and is therefore 
described by an angle < < 27r, and its conjugate angular momentum 
Pi, with i = 1,. . . ,N. The constant J > is the interaction strength. The 
aXY Hamiltonian has been extensively studied in the last few years. It has 
been shown that the long-range cases {a < d) have the same thermodynamic 
properties as the a = case (HMF), in all the energy range [11,12]. If one 
looks only at the high energy regimes, then for all a the equilibrium properties 
are those of a gas of almost free rotators. However, time scales do depend on a. 
This is the case of LLEs, which have been studied in detail, both numerically 
[10] and analytically. Scaling laws for high energies were obtained by Firpo and 
Ruffo [13], using the geometric method [2,3] and also by the present authors, 
using a random-matrix approach [14]. 

In order to apply the stochastic approach, in its mean-field second-order- 
perturbative version, one has to calculate the average, variance, and corre- 
lation function of the Hessian V(t), i.e., Eqs. (10) and (11). For high energies 
these ingredients can be easily estimated, using simple approximations. Any- 
way, we will verify that the theoretical estimates, coarse as they are, coincide 
with the exact time averages calculated numerically. 

The elements of the matrix V can be readily obtained: 
J cos{9i - 9j) ^ 
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Then, inserting this result in the definition (10) of n we arrive ai ji/J — 1 — 
2(V) / {NJ). The average potential energy (V) can be calculated as a canonical 
average following the (now standard) saddle-point calculations of Refs. [12,13]. 
Skipping the details, one arrives at 



where (3, the inverse of the temperature, is related to the energy per particle 
e through (3~^ = 2e — J, in the disordered regimes. This expression for fi 
coincides with the result obtained by Firpo and Ruffo [13], except for a factor 
of two. 

The calculation of proceeds in the same way as in the case a = (see 
[7]), but here we make the extreme approximation of free rotators. That is, 
we assume 



(cos(^^,-^^,))^0, 

{cos{9i - 9j) cos{9i - 9k)) ~ , (15) 

where it is understood that i ^ j and j ^ k. Then the final approximate 
expression for o"^ is 

N{a,d) ^ ' 



At high temperatures, the relative importance of the interactions decreases 
with increasing energy and N ^ and the dynamics is dominated by the kinetic 
part of the Hamiltonian. The picture is that of particles rotating almost freely 
during times which are long as compared to the mean rotation period. For 
any value of a we can assume, as a first approximation, that the values of 
{6'fe} and {pfc} are independent random variables with uniform and Maxwell 
distributions, respectively. The decorrelation mechanism is just single-particle 
phase mixing. Within this approximation the correlation function is Gaussian 
and the correlation time is of the order of the mean period of rotation, i.e., [7] 




Gathering the results of previous sections, one can construct the 3x3 matrix 
A, and extract the LLE from the eigenvalue of A with the largest real part. 
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However, we have verified that for the cases considered here, the LLE is very 
well approximated by the asymptotic expression [7] 



A 



1/3 



(18) 



with cr^ and given by Eqs. (16) and (17). The absence of in the leading- 
order expression for A says that in the disordered phases of the aXY , fluctua- 
tions are much larger than the average, i.e. a ^ fi, and dominate the tangent 
dynamics. We conclude this section by noting that Eq. (18) leads to the same 
scaling laws obtained in Refs. [13,14]. 

4 Numerical studies 

First of all we verified that the analytical calculations of /i, cr^, and Tc agree 
with the corresponding dynamical averages. This test was perfomed for the 
case d = 1. In Fig. 1 we show numerical values for fi and for various values 
of a. Numerical integration of the equations of motion was performed using a 
fourth order symplectic algorithm [15], with equilibrium initial conditions, i.e., 
angles uniformly distributed in [0, 2n) and Gaussian distribution of momenta. 
Forces were calculated using a Fast Fourier Transform algorithm. One sees 
that the numerical data can be well described by the approximate expressions 
(14) and (16). 



10° F 




10.0 



Fig. 1. Numerical (symbols) and analytical (lines) values of fi, cj^, and Tc for the 
aXY model as functions of the interaction range a {d = 1). Energy is fixed at 
£ = 5.0, J = 1 and system size is N = 128 (circles) or N = 512 (squares). Error 
bars (not shown) are of the order of the symbol size. 

In order to discuss the correlation times, let us now look at the correlation 
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functions of Fig. 2. Theory and simulations agree almost perfectly in the cen- 
tral part of the distributions. However, there are long tails that grow with 
increasing a. We recall that the quantities t^^^ are moments of the correlation 
function /(r), and as such, very sensitive to the precise shape of the tails. So, 
as a grows, we lose confidence on the Gaussian estimates for t^''\ 
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Fig. 2. Correlation functions for N = 256, e = 5.0 and J = 1. Several values of a are 
considered (d = 1). We have averaged over 100 initial conditions. The dotted line is 
our theoretical prediction (a Gaussian) . The inset shows that tail fluctuations grow 
with a. 

For comparing the theory with simulations, we will use LLE data existing in 
the literature [10,16]. Figure 3 shows LLEs for different values of a as a func- 
tion of the system size and £ = 5.0. We have considered large particle numbers 
(A^ > 50) and an energy well above the transition to ensure (i) the validity 
of the approximations invoked in calculating the averages of Sect. 3, and also 
(ii) to guarantee that we are in a disordered, quasi-bahistic regime. Full lines 
correspond to the theoretical LLEs obtained by using the asymptotic expres- 
sion of Eq. (18). Figure 3 shows a satisfactory agreement between theory and 
simulations. Differences in the definition of N were absorbed in J, otherwise 
set equal to one. 

5 Summciry and concluding remarks 

In a previous paper we had seen that the stochastic recipe works satisfactorily 
in the quasi-ballistic regimes of the HMF (a = 0). Now we have verified that 
the theory still works well when the interaction range is not infinite. Even for 

short-range interactions we have found a good agreement between theory and 
simulations. This suggests that the crude approximations that allowed us to 
estimate the LLE may be reasonable in these cases. Let us summarize the 
approximations we used. 
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Fig. 3. Largest Lyapunov exponent of the aXY Hamiltonians as a function of system 
size N, interaction range a and lattice dimension d. Energy is fixed at £ = 5.0. 
Symbols correspond to numerical data in the literature for d = 1 [10] and d = 2,3 
[16]. Lines are our theoretical predictions, labeled by the ratio a/d. 

First of all we invoked a weak intermittency approximation to exchange log- 
arithm and average. We tested this approximation in the a = case [7] and 
verified that it does not introduce an important error (about or less than 10%, 
for the cases considered). Even though we have not made numerical tests for 
q; > 0, we suspect that intermittency does not increase with a, one reason be- 
ing that the system becomes more chaotic with increasing a (the LLE grows 
with a). 

The second approximation, and most important one, is the truncation of the 
cumulant expansion at the second order. By doing so, one introduces a relative 

error of the order of the Kubo number aTc. This perturbative parameter grows 
with a/d but saturates at a finite value as « ^ oo (as illustrated in Fig. 1 for 
d = 1), so there are no indications that the cumulant expansion may fail for 
short-range. 

The "mean-field" diagonalization is exact in the infinite-range case [5] and was 
expected to be a good approximation for long-range interactions {a < d). For 
short-range forces, this approximation represents a crude truncation of the 
basis for diagonalizing A. Note however that we are dealing with a Hermitian 
problem, i.e., finding the largest eigenvalue of {$,^'^){t) [Eq. (7)], then, if A is 
calculated accurately, truncation of the basis may still produce a good lower 
bound to the LLE. This could be the explanation for the agreement between 
theory and simulations also ior a > d (see Fig. 3). However, this accord should 
be further explored and tested for values of a larger than those here considered. 

We have limited our presentation to equilibrium situations. If the system does 
not explore the full phase space, as is the case of the quasi-stationary regimes 
of the aXY [17], the stochastic method may still be applied, but then one has 
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to replace microcanonical averages by averages over the (a priori unknown) 
regions of phase space that are effectively visited by the system. 

We are grateful to A. Campa, A. Giansanti, D. Moroni and C. TsaUis for com- 
municating their numerical data. We acknowledge Brazilian Agencies CNPq, 
PAPER J and PRONEX for financial support. 
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